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1. Introduction 

One of the most interesting open fields in quantum magnetism is tlie study of 
frustration [1] . Frustrated models are described by Hamiltonians with competing local 
interactions such that the ground state cannot minimize their energy simultaneously. 
Frustration can appear due to the geometry of a problem, the Ising model on a 
triangular lattice being a paradigmatic case, or due to the coexistence of ferro- 
and antiferromagnetic interactions, as in spin glasses. Frustrated models typically 
have highly degenerate ground states, which can become ordered by increasing the 
temperature or by quantum fluctuations — i. e. "order by disorder" . Theoretical and 
numerical problems, such as the large dimensionality or the sign problem in Monte 
Carlo simulations make it very difficult to study frustrated Hamiltonians. 

Quantum simulation has been put forward as a tool to probe the physics of a 
variety of many-body systems. In particular, schemes with cold atoms in optical 
lattices have been suggested to simulate arbitrary spin models [21 [Sj IH El [6l [7] and 
some interesting frustrated Hamiltonians [51[S] . Compared to the alternative of looking 
for or inducing low-dimensional behavior on existing magnetic materials [ 101 lll| and 
in organic conductors [12, J3J , the cold atoms and molecules offer greater flexibility in 
terms of variable geometry and interaction strength. However, the majority of current 
proposals are perturbative and their effective interactions are rather weak. This makes 
their experimental realization challenging, as it requires very low temperatures. A 
solution suggested to solve the problem of weak interactions is to replace the atoms 
with polar molecules [TT. 

In this paper we follow a different approach. First of all, spin up and spin down 
states are identified with particle and hole configurations. This gives naturally an XX 
coupling which is proportional to the tunneling amplitude J15j and, thus, it is strong 
and less experimentally demanding. By using the internal state of the atoms to encode 
a second virtual dimension and applying a spatially dependent Raman coupling, we 
are able to induce in-plane frustration. The result is a new frustrated model, which 
contrasts with related literature on spin ladders, where the Hamiltonian either use 
antiferromagnetic isotropic [16l [171 [HI [HI O [21] or XXZ |22] interactions, or use 
purely XY interactions but require a different lattice geometry '^, or both. More 
important, the problem we consider here is a full Bose-Hubbard Hamiltonian, and 
the main result is that the transport properties of the lattice change dramatically, 
causing the breakdown of superfluidity even for large densities in which the analogy 
to spin models no longer applies. In particular, in the most interesting case the Mott 
insulator is replaced by a new dimerized phase made, where neighboring pairs of sites 
have maximal coherence and there is a fast decay of coherence between pairs of sites 
as a function of distance. This new phase is reminiscent of a Bose glass [24j, with 
the important difference that it has been induced by frustration and not by disorder. 
Our statements are supported by a variety of analytical, variational and numerical 
solutions. 

The paper is organized as follows. In Sect. [2] we introduce the two Hubbard 
models that we are going to use, explain how they are implemented with cold atoms 
in optical lattices and briefly justify the use of the word "frustration". In Sect. [3] we 
focus on a model with diagonal interactions across each square plaquette. We first 
discuss the different phases as obtained from DMRG calculations, including qualitative 
arguments about why these phases are expected. Afterwards we explain in more detail 
how fragmentation happens and what happens in the limit of strong interactions, and 



Fragmented superfluid due to frustration of cold atoms in optical lattices 3 

what are the correlation properties of all available phases. In Sect. 3] we briefly present 
a model with frustrated rectangular lattice. In Sect.[5]we study in great detail all issues 
regarding the implementation of our ideas in current experiments with optical lattices. 
First we do a microscopic derivation of the parameters in the Hubbard Hamiltonians, 
including all possible sources of error. The main conclusion is that the requirements for 
studying these frustrated models appear to be within reach of current experiments and 
that there are no side-effects from introducing a Raman coupling in the experiment. 
We explain how the different phases studied in this paper could be detected, either 
from time of flight measurements or more sophisticated correlation measurements. 
Finally, we discuss possible sources of imperfection, such as temperature or a residual 
harmonic conflnement. The last section (Sect. [6]) contains a brief summary of our 
main results and possible implications and connections to other works. 

2. The model 

2.1. Bose-Hubhard m,odel 

The system that we are studying is that of cold atoms confined by an off-resonance 
optical lattice that forms a ID bosonic lattice gas [15]. As in recent experiments [25], 
the lattice traps atoms in two internal states. If the conflnement is strong, the effective 
model will be a Hubbard Hamiltonian [26] 

H ^ -J '^ a\aO-3a + ^ -^n^^Uic,' + Hb., (1) 

{ij),a iaa' 

where J is the hopping amplitude and we assume that the on-site atomic interactions 
are both repulsive and symmetric, Uaa' = U > 0. Note that in this work, the internal 
state of the atoms are denoted by Greek letters, a =t, i; the site indices are denoted 
by Roman characters, i, j, and go from 1 to L which is the length of the lattice; flnally 
N := N^ + Ni, represents the total number of particles. 

The new ingredient in our Hamiltonian are two counter-propagating laser beams 
that induce Raman transitions between the two atomic states. By adjusting the phase, 
the polarization and the alignment of these beams, we can ensure that the effective 
Rabi frequency, ^{x), forms a lattice with twice the period of the conflnement [Fig.[T]. 
We can interpret our one-dimensional system as having two ID lattices, one for each 
internal state of the atom. Both chains are coupled forming a ladder thanks to both 
the interaction, U, and to Hn, a Raman term with the spatially dependent Rabi 
frequency, Q{x). 

We shall consider two specific configurations of the Raman coupling. When its 
maxima and minima coincide with those of the confining lattice [Fig. [TJi] , the result 
will be on-site Rabi oscillations 

i/7? = 5](-l)V'(a,\aa+H.c.). (2) 

i 

Alternatively, when they are displaced by half a period we will have diagonal 
interactions [Fig. [1^] 

Hr = 5](-l)V(at^i^a.i + al+^a^T + H.c). (3) 
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Figure 1. (a) Raman-assisted tunneling. A second optical lattice (dash-dot 
line) connects neighboring Wannier functions of atoms in different states with 
amplitude J'. On the other hand, normal hopping moves atoms while preserving 
the state, J. (b) The hard-core bosons limit resembles a spin ladder with 
frustrating diagonal interactions, (c) The change of variables in Eq. (O replaces 
the frustrating terms with weak and strong tunnelling, (d-f) Similar drawings for 
the Villain model in Eq. ^. 



2.2. Origin of frustration 

The frustration induced by Hji becomes evident in the hniit of hard-core bosons, 
af^ = |15| . Identifying the bosonic operators, a^j and a^i, with the spin operators 
(j~-^ and (J~2, the Hamiltonian ([1]) becomes an XX model on a spin ladder [Fig. [TjD,e]. 
Thus, while the hopping J translates into a ferromagnetic XX interaction along the 
legs, the Raman coupling H^i is a transverse XX interaction which alternates from 
ferro- to antiferromagnetic. 

Take for instance the diagonal interactions ([3|). The associated spin model is 



H = 



E 



E 

--l.2j=l...N 

E 

= 1...7V 






ij^ij + 1 



V V 



l) + 



= 1...JV 



(-l)V'(af,ar,+i 



^■'^Ij+i) 



y V 



(4) 



The frustration arises from the competition of positive and negative contribution 
terms. For instance, looking at the sites (l,j),(l,j + l),(2,j + 1) and (l,j + 2): 
the bonds on this triangle have three ferro- and one antiferromagnetic terms, and 
there is no way in which the spins can be aligned so as to minimize the energy of all 
terms [see also Fig. [TJ)]. 

While the previous analogy is somewhat pleasing, one may wonder whether having 
more than one particle per site or working with a system which is not properly 2D 
will give rise to a trivial Physics. Regarding the first point, when we have more than 
one atom per site our model becomes equivalent to an array of Josephson junctions 
[9] where the frustration is still present and lays on the choice of phase for each 
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individual site. Regarding the second point, some of the best studied models in 
frustrated quantum magnetism, such as the Majumdar- Ghosh [571 HH] and the zig- 
zag chain [5^ are also quasi-lD models. As we will show in the following sections, 
the competition between the different terms in Eq. ([1]) does indeed produce new and 
interesting effects with the advantage of being analytically tractable. 

3. Diagonal interactions 

In contrast to fully two-dimensional systems, in our quasi- ID models can be greatly 
simplified by rewriting the frustrated coupling using dressed states 

Cz±-=—^{aiT±a,i). (5) 

We will begin by studying the Raman interaction ^. With the previous change of 
variables, our Hamiltonian becomes 

^ = - E [^ + ^(-l)V'](4c,+i,-HH.c.) + |n2, (6) 

2;cr— ± 

where n^ :— cl^Ci+ + c\_Ci- is the total population of a single site. Eq. ([6]) models 
a ladder whose legs are made of junctions with strong, J+ = J + J' , and weak 
hoppings, J_ = I J — J'\, as shown in Fig. [T]:. Notice that the model is symmetric 
under the exchange of J and J'. In the following we will often exchange descriptions 
between J+, J_ and J, J' , and we will assume, without loss of generality, that all these 
parameters are not negative. 

3.1. Phase diagram. 

In this section we present zero temperature results for Hamiltonian ([5]). We have 
computed the ground states and first excitations of this model using the Density 
Matrix Renormalization Group (DMRG) method in the Matrix Product States (MPS) 
formalism [29l [30]. We have conducted accurate simulations for L — 8, 16, 32, 50 and 
70 sites, using a cutoff of 4 particles per spin. These values give us a local dimension 
of the Hilbert space of 25 states, which is very big and can only be handled with 
state-of-the art optimizations, such as working in sectors with well defined number of 
particles or angular momenta. Due to the large size of the local Hilbert space, we have 
typically worked with MPS of size 100 and checked convergence for different points 
with up to Z? = 20C|j. In our simulations we have computed the ground state and 
lowest energy states for N/L = 1/2, 1, 5/4, 3/2, 7/4, 2 particles, as well as the energy 
to add a particle, ^p :— E'a'+i — E'at, or to create a hole, ^h '■= Epf — -Eat-i, and 
verified that these values did not change with a larger cutoff. 

These zero temperature simulations reveal a very simple picture that we 
summarize here. As illustrated in Fig. [TJi, when we increase U for fixed J_|_ and J_, 
the system first experiences a phase transition from a superfluid to a fragmented or 
dimerized phase in which coherence is maximal between pairs of neighboring sites. As 
the interaction is further increased these fragments evolve smoothly (i. e. a crossover) 

I Note that due to our way of computing the two lowest excited states, which uses diflferent MPS 
for the ground state and excitations. These simulations are comparable to DMRG calculations with 
D = 300 and D = 600 states, respectively, because, in the DMRG parlance, the effective basis states 
are optimized independently for each excited state. 
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Figure 2. (left) Phase diagram for the diagonal interactions JHJ computed 
using DMRG. For J' = 0.9J, the gray areas denote the incompressible Mott 
lobes in the space of chemical potential (fi/U) vs. hopping amplitude (J/U). 
The in between regions are in superfluid phase. For J' = J, the phase space 
fills with incompressible regions, each one having integer or half-integer filling, 
N/L = 1,3/2,2,.... These regions are delimited by the dash-dot lines, (right) 
Schema of the different phases as a function of the interaction U. 



to a Mott insulator. While the first phase is gapless and its excitations are phonons, 
the crossover contains a gapped phase whose lowest excitations are localized spin flips. 
In the following sections we will study in detail the properties of the different 
coupling regimes, numerically as well as analytically where possible. Nevertheless it 
is possible to obtain a qualitatively picture of the three regimes and the associated 
ground states by variational wavefunctions. First of all, for weak interactions, U <^ J±, 
the ground state is a uniform superfluid that spread over the lattice 



IV'Sf) oc 



1 



^Nl 



E 



{aal^ + (3al 



il' 



N 



(7) 



with the usual rotational symmetry on the (a, (3) space. 

For stronger interactions, J_ < C//4 <C J+, the energy of IV'Sf) is larger than 
a state made of fragmented condensates which reside on the junctions with large 
hopping, J+ [Fig.[TJ:], 

t 



t 

2k+l, 



-'2k+2, 



r- 



vac) 



(8) 



l^frag) OC n(4/c.+ +4^+1.+ )"+ X (ci 
fc 

Here, for integer filling, N/L = 1,2,3,..., the ground state has is one with 
spontaneously broken rotational symmetry, a ferromagnetic state with either n+ = or 
n_ = 0. This effect arises from the contributions to the energy of both the interaction 
U and the hopping J_, as it will be explained in Sect. 13.21 

Another symmetry break happens when the interactions become dominant, 
U ^ J±. Indeed, in the limit of strongly interacting bosons and N = L, all ground 
states can be written in the form 



|^Mott)(xn(4+)"'+(cL)"'ivac) 



(9) 
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Figure 3. (a,b) Finite size scaling of the energy gaps for spin and density 
excitations, AEi and AE2, for a particular case J' /U = 0.9, J/U = 2.8. We 
fit both curves to a polynomial in 1/L for L = 4, 8, 16, 32, 50 and 70 and take 
the limit 1/L -^ 0. With this technique we plot in (d,e) the energy gaps 
vs. hopping amplitude, for J' /J = 1.0,0.9,0.8,0.6 (top to bottom). With 
vertical lines we mark the points where the breakdown of superfluidity is roughly 
expected to take place, C/iij ~ 3.84|J' — J\. (c) Mean correlation length, 

f ~ Y.i,s W'^U^+s)/ Y.i,s EK^ «»+«)' fo'- J' /J = 0-6, 0.8, 0.9, 1.0. 



The quantum fluctuations generated by the hopping select, among aU disordered 
ground states, one particular ferromagnetic state with either n_ = or n^ = 0. 
This case is further discussed in Sect. 13.31 

In the space of parameters, J/U, vs. chemical potential, /i, we can draw a phase 
diagram similar to that of the Mott-insulator transition 24 . As Fig. [2] shows, for 
J' — 0.9 J, the curves of particle and hole excitations, fip and ^h, create lobes (gray 
areas) containing fragmented phases. Crossing these lines amounts to a quantum 
phase transition to a superfluid, on the exterior of the lobes. Qualitatively, these 
lobes are similar to those of the Bose- Hubbard model [24,, but about J/\J — J'\ times 
larger, so that for ,]' = ,] the lobes become inflnitely large and all ground states are 
fragmented. In this limit of J' = J and for integer and half-integer flUing, that is 
N / L = 1, 3/2, 2, . . ., we have on phase space infinitely many incompressible "stripes" 
whose borders, ^p and ^/i, can be computed exactly [See Sect. 13.21 and Eq. (H^)]. As 
explained above, within each stripe the state of the atoms changes smoothly from a 
Mott insulator (small J+) to fragmented or dimerized phase where atoms delocalized 
between pairs of neighboring sites (large J+). Interestingly, in the macrocanonical 
ensemble, all phases with other filling factors collapse to the lines between these 
regions. The large degeneracies of these other fillings has made it impossible for 
us to study systematically the nature of the those states and of their excitations. 

To gain further insight on the properties of the ground state, let us focus on the 
unit filling case. The transition from the Mott phase to the superfluid is accompanied 
by the change of several excitation gaps. One is related to the energy required to 
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add or remove a particle, the other one is the energy gap required to create a density 
perturbation for fixed number of particles. Both gaps close at the same point when 
we jump from the fragmented phase to the uniform superfluid, signaling a phase 
transition from an incompressible to a compressible phase and also the establishment 
of algebraically decaying correlations [Fig. [3] . This transition occurs around the critical 
value of the interaction for a single-component ID Bose-Hubbard model [3T | [32 l [33 l [34] 
with hopping J_, that is f/ ~ Uid '■= 3.84J_ (See Sect. 13. 4|) . A third energy gap is 
present for excitations that take a particle out of one of the coherent fragments into 
a different spin state, while preserving the total number of particles. This gap freezes 
the dynamics of the spin of the atoms. As shown in Fig. [3^-b, it is maximal in the 
fragmented phase, with a maximum value of order C//4. For strong interactions it 
grows as J'^/U and in the superfluid phase it should close. However being U ^ Jpm a 
region gapless with respect to density excitations, the DMRG simulations are probably 
not accurate enough to reflect this fact. 

3.2. Fragm,ented state 

In this subsection we will comment on the configurations that appear for U G [J_, J+). 
For simplicity we will begin with the case of balanced hopping and Raman coupling 
and no interaction, [/ = J_ =0. Then it is trivial to see that all ground states will be 
of the form ([5]) with atoms delocalized on pairs of sites. In particular, for unit filling 
N — L we have three degenerate states: a uniform state with all junctions populated 
n± — 1 [Fig. \^] and two states with broken translational symmetry, n_ = 2 and 
71+ — 2 [Fig. [4^]. This degeneracy persists for small U, but it is broken as soon as we 
set a weak imbalance in the hoppings. For any < ./_ ^ U, terms proportional to J_ 
connects these states ofF-resonantly to excited states like the ones shown in Fig. |3j:-d. 
It follows that to second order in perturbation theory the degeneracy of these states 
is broken and the energy of the uniform state, n_|_ ~ n^, is shifted by an amount 
AE — —8{J'^/U) X L smaller than the ferromagnetic states with broken symmetry, 
AE — — 24(J^/C/) X L, which become the true ground state. 

Once we have established that the system spontaneously selects either the c+ or 
c_ order, we can compute the ground state for any value of U while keeping J_ — 0. 
Such state will be a product of ground states of the two- well problem, 

i/2 n 

1^) = J] |^(«)) = n E 4,?(4.J''"(4+i.)''"'"|vac). (10) 

i— 1 i—l m—0 

Here fi = N/L = 1 is the density, (f>m is the wavefunction of the double-well problem 
with n atoms and cr = ± is the selected polarization of the atom. 

To delimit the incompressible regions in phase space we only need to compute 
the energies of a state with N = L,L + 1 and L — \ particles, which are variations of 
the one in Eq. (|10p . Since the ground state energies of a double well with 1, 2 and 3 
particles are, respectively, 

ei(J+,[/)= -2J+, (11) 

e2(J+,f/) = |- 
e3(J+,C/) = 2C/- 
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Figure 4. Two possible ground states for tiie fragmented phase: (a) witii and (b) 
witiiout spontaneous symmetry breaking. A nonzero value of J_ couples the state 
(a) and (b) to the excitations (c) and (d), respectively. This coupling makes the 
state (a) preferred for any nonzero J_ in the fragmented regime J- < U < J+ . (e) 
A particle hopping through the lattice sees a kind of Mach-Zender interferometer. 
The antiferromagnetic hopping (dashed) gives rise to destructive interference on 
all paths longer than 2 sites. 



the boundaries of the Mott region with n = 1 particle per site are given by 

Aip(J' = J)=e3(2J,C/)~e2(2J,C/), (12) 

Aih(J' = J)=e2(2J,C/)-ei(2J,t/). 

These boundaries are the dash-dot hnes we plotted in Fig. [2l A simple inspection of 
the previous formulas shows that these curves never cross and become parallel in the 
limit oi J+/U > 1. 

The previous exact results, which are valid for J_ — 0, can be extended to 
J_ <ti U, J4- by means of perturbation theory. In that case there will be an additional 
coupling between fragments up to a distance of the order of the correlation length. 
The wavefunction can be approximated by linear combinations of the terms ()10|) with 
different occupations of the double wells. The corrections are of order 0{J^/U) to 
the wavefunction, and of order 0{J^/U) to the energy, but the dominant term is 
always the one with equally populated double wells. This will be further explained in 
Sect. 13. 4[ where we analyze the location of the Mott insulator to superfluid transition. 
That this picture is applicable to a great degree is appreciated in Fig. [5l where we 
plot the correlations between different sites. For J' = 0.9, sites connected by J+ 
are strongly correlated, while the correlations decrease significantly between different 
pairs of sites. 



3.3. Hard-core bosons 

To understand the behavior of the bosons in the strongly interacting limit, C/ 3> J± it 
is useful first to see what happens if there are no frustrating terms, J' = 0. We know 
that on the one hand, if we are below half-filling, N < L, the gas is a Tomonaga- 
Luttinger liquid with spin-charge separation [35]. In the limit of C/ — > 00 the spin 
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J/U=0.6 




Figure 5. Correlations between diflFerent sites, Cij := (cj, Ci-|_), for J/U = 0.6 
and L = 32 sites. In solid line we plot cig j for different values of J' /U. For 
the strongly coupled sites 15 and 16, the correlations are large. For J ~ J' the 
correlations are dramatically reduced between sites that belong to different double 
wells in the effective superlattice. Note that the correlation function has period 



two, Ci 



■ Ci_|_2,j-|-2, due to the superlattice structure. 



degrees of freedom become degenerate and the whole system behaves hke free fermions 
or a Tonks gas. On the other hand, for half filling, N = L, the sample becomes a Mott 
for any strong interaction U ^ J. The atoms cannot move but their spin degrees of 
freedom interact with a ferromagnetic Heisenberg interaction of strength t ~ —,P/U 
and a gapless spectrum made of spin waves. 

The diagonal coupling J' changes this landscape slightly. For half filling, N — L, 
we still have a Mott phase, with one atom per site. Connecting with the previous 
section, as we increase U the wavefunction of the two-well fragment PH)) becomes 



closer and closer to 



a(«) 



ii/2- 



This is a smooth evolution which is why we 



speak of a crossover. Furthermore, compared to the unfrustrated case J' = 0, the spin 
rotational symmetry is now broken, favoring the X direction or c± states in the second 
quantization language. The effective spin model is, up to constants, a ferromagnetic 
one 



1 



^snin 



^spm 



^r 



hj'^aza 



k"k+l 



J+J-ffkO^k+l] 



(13) 



and spin excitations are prevented by an energy gap A = 2J'^/C7 |36j . This gap is 
evident also in Fig. [3^, in the region of small hopping, J ^ U. 

For smaller fillings, i. e. N < L, the previous spin model is no longer valid 
because the chain has holes. Nevertheless one can easily compute things if J_ = 0. 
The ground state manifold is composed of variations of the fragmented wavefunction 
(jlOp . In particular, the set of available configurations is shown in Fig. [6l where we 
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Figure 6. (a) A possible state of the ladder, where the circles denote the location 
of a particle, (b) The same state mapped to a ID problem for J_ = 0. Allowed 
(b)-(c) and forbidden (d) states in the hard-core bosons limit {U — > oo). 



show three kinds of bonds. Those with zero or two particle, which are frozen, and 
bonds with one particle on the left or right sites, which can be identified with spins. 
These effective spins form ferromagnetic blocks which are separated by the frozen 
sites. If the block is surrounded by at most one doubly occupied bond, then it can 
only adopt one state and it has zero energy [Fig. [6jl] • If the block, on the other hand 
is surrounded by empty blocks it may host at most one domain wall [Fig. [61;] with 
some small negative energy due to the motion of the domain wall (for instance, the 
states in Fig. [Sb and Fig. [St are connected by hopping). The effective Hamiltonian 
for a block of size M is 



M 



HDw^-jy(\m+l){m\ + \m){m + l\), (14) 



where m denotes the number of particles to the left. This gives an energy spectrum 
€k{M) = — 2 Jcos[fc7r/(Af + 2)], for fc = 1 . . . M+1. Since efe(M) is a decreasing, convex 
function of M, the lowest energy state is achieved by making L — N + 1 blocks with 
about N/{L — N + 1) atoms each. The excitations are still gapped with A ex J, and 
they can be of two types: excitations of the domain wall, and merging of two blocks. 
If A^ < L/2 the latter are more important and of order A = 2ei(l) — ei(2) ~ 0.5858 J. 
On the other hand, if A^ > L/2 the excitation energy of a domain wall can be pretty 
small, 0{l/A'P). 

3.4- Scaling of the Mott transition 

In this subsection we attempt to analyze how the superfluid to insulator transition 
changes due to frustration. We will do it in two different perturbative methods. The 
first one is based on a strong coupling expansion around J, J' — 0, while the second 
one is based on a strong coupling expansion around J' — J. Remarkably, this second 
method, which is more accurate than the trivial strong coupling expansion 33J, may 
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Figure 7. Critical value of U/J at which the Mott insulator to superfluid 
transition should take place. The blue and violet curves are based on the strong 
coupling expansion of third and second order in {J/U) developed in Ref. 1331 . The 
solid line is based on a different strong coupling expansion of order 0{\J — J'\'^/U), 
developed below. The dotted line is the critical value Uc = 3.84J computed with 
DMRG for unfrustrated models 32 . Note the agreement between the two latter 
methods for J' = 0. 



be contemplated as a real space renormalization technique in which the ahcrnating 
Hubbard model is replaced by an effective one with a renormalized hopping and half 
the sites. 

Our first attempt starts with the estimates of the critical value of J/U given in 
Ref. [33]. Applying perturbation theory up to 0{J^/U^), it is possible estimate the 
change in the energy induced by a small amount of hopping. The resulting energies 
can then be used to compute the particle and hole excitation gaps, and delimit the 
borders of the incompressible regions in the hopping vs. chemical potential parameter 
space, as in Fig. [21 To apply the theory in Ref. [33 we need the hopping matri4§l 

— tij = —JiSi+ij — Ji-i5i^ij, (15) 



the smallest eigenvalue of {-—tij) 



K 



= -max{|J+ - J_|,|J+ + J_|} = -2max{J, J'}, 



(16) 



the wavefunction of the single-particle ground state, fi = 1/v L, and finally the sums 



E4 



fi = J 



Ji 



ij 



Jitijjj — J+ + J- 



(17) 



As shown in Fig. [3 the transition does indeed shift to larger values of J/U as we 
switch on the diagonal coupling J'. However, the prediction is not very consistent 
between second and third order expansions, and it even fails to provide the right value 
for J' ~ 0, probably due to the asymptotic nature of these expansions. 

Our second attempt focuses on the regime of strong frustration or values of J' 
close to J. We have seen that for J' — J the incompressible regions fill the entire phase 

§ We assume periodic boundary conditions. 
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space. With very simple arguments it is possible also to compute how they shrink for 
nonzero \J' — J\ = e -^ U. The reasoning begins by noticing that our system is a 
superlattice, where the weak hopping terms can be treated perturbatively. However, 
instead of using the Mott states as unperturbed states one must use the exact ones for 
J+ > 0. Basically, we have an unperturbed state with energy E = Le2/2 and several 
unperturbed manifolds above it, consisting on holes and excitations in the superlattice, 
with energy gaps of order 63 — 62, 61 — 62. We can compute the corrected fip and Hh 
as a function of J_ . The equation for the critical point then looks as follows 

Mp - Mh = 63(C/) + 6i([/) - 262(C/) - J-ciU) + OiJl) = 0, (18) 

where the energies and the matrix elements of the hopping Hamiltonian, c(U), are 
functions of U. We have solved numerically the highly nonlinear equation P^ for U 
as a function of J_ and J_|-. The results are plotted in Fig. \7\ As seen there, this 
second method interpolates properly between the critical value of C/ = 3.84 obtained 
by DMRG for the unfrustrated model [32] , and the critical value U = computed for 
J — J'. Furthermore, the whole line deviates very little from Uc — Uc{J' — 0)| J — J'|. 

3. 5. Correlators 

As it will be useful to identify the different phases experimentally, we have computed 
the single-body and two-body correlation functions in the three regimes. The most 
interesting values are the single-body correlators, {ajaj). In the superfluid limit, due 
to the one-dimensional character of our problem, we expect that correlations decay 
algebraically. In the fragmented regime, we have two possible ground states, depending 
on whether the fragments sit on even or odd junctions, which is equivalent to specify 
whether the fragments are made of c-|_ or c_ particles. In the first case we will have 

(4,, + C2», + ) = (4,, + C2.+l, + > = (4, + C2*, + ) = 1/2, (19) 

with all other correlators being zero, and in the c_ case we will have the same 
correlation pattern but displaced by one lattice site. Finally, in the Mott insulator 
regime we have no correlations and (c|^Cj±)frag ~ ^i-j- 

In experiments what is measured is the time of flight images. These pictures are 
related to the momentum distribution of the sample or the Fourier transform of the 
single-particle correlations 

L 

n, (X ^ e'^.'C^-jVL^^t^^.^ ge[-L/2,L/2]. (20) 

k,l=l 

This function is peaked at q = and decays towards the sides of the Brillouin zone. 
The half-width of the peak is related to the correlation length of the sample. In 
the superfluid region it will be proportional to the size of the lattice and limited by 
finite temperature effects |37i . As soon as we cross the phase transition towards the 
fragmented phase, though, the correlation length decays abruptly to ^ = 2 sites and 
the previous function has a constant width 

n, - 1 + cos{2TTq/L) (21) 

From there on we expect a smooth crossover towards riq = \ which is the Mott- 
Insulator regime. All this phenomenology is evident in Fig. [512 
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4. Villain model 

The other Umiting case in our setup is that of vertical frustrating interactions, a 
configuration resembling the odd model or Villain lattice [1] . Once more it is convenient 
to change basis as in Eq. ([5]), so that the frustrating terms become alternating energy 
shifts 



^-E 



- J(4c,+i, + H.c.) - (-l)VJ'cLc,, + -n^. 



(22) 



i;cr— ± 

The effect of the Raman coupling is now equivalent to a superlattice. This additional 
potential splits the Wannier band into two effective bands separated by a gap J', 
each band having a reduced width AJ := (-\/4J^ + J'^ — J')/2. The superlattice 
localizes the atoms on alternating rungs, forming an antiferromagnetic configuration. 
Instead of exhibiting fragmentation, the system goes straight from a superfluid to a 
Mott-insulator, but now the transition happens for weaker interactions, U ~ A J. 

5. Experimental realization 

There are several issues that one has to consider when implementing our Hubbard 
models using neutral atoms. The first one is how to relate U, J and J' with 
experimental parameters such as the intensity of the Raman laser and the strength 
of the optical lattice. Additionally one has to consider how to prepare the ground 
state and how cold the sample should be. Finally one has to think about procedures 
to detect the different phases. We will address all of these questions in the following 
sections. 

5.1. Microscopic theory of Hubbard model 

In this subsection we show how to derive the constants U, J and J' using a 
band structure calculation that takes into account the Raman coupling, n{x) — 
flo cos{kx + 0), and the confining lattice, 

V{x) = Vo^sm{kxf + Vo^[s[n{kyf +sm{kzf]. (23) 

Our work generalizes the ideas from Ref. [26] , but now our microscopic Hamiltonian 



contains a potential term 



-|lv2 + M/..,(x) 
Zm 



Ua 



V^.'(x) (24) 



W,,, = V{^)S,,, + n{^)S,,, + Y, ^l^c(x)|2 (25) 

a=T,i 

that includes not only the trapping potential, but the coupling and the interaction 
between atomic states. 

Following |26j , we perform a tight- binding approximation and replace the bosonic 
operator by an expansion of the form 

V'(x) ~ ^ aijkaWx{x - bi)wy{y -bj)wz{z - bk), (26) 

i,j.k,(7 

where b = ir/k is the period of the lattice and Wx,y,z are the Wannier function 
corresponding to the lowest Bloch band of the lattice, along the X, Y and Z directions. 
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respectively. By substituting this expansion into the full Hamiltonian ([M| and 
retaining the most important terms, one arrives to a Bose-Hubbard model H]). 

To compute the parameters of our model, we will assume that the confinement 
along the transverse directions, Vo_l, is big enough, so that hopping is negligible 
and our system is made of multiple, disconnected ID lattices. However, even with 
this assumption there are still too many "knobs" to tune. To simplify the problem 
it is customary to adimensionalize all magnitudes using the energy of the lattice, 
Er = fi fc^/2?Ti, and its period, 6, as units. We replace the trapping potential 
with V ^ Vb cos(2fca;), make the change of variables 2kx — > 2ttx and rescale the 
wavefunction accordingly, w -^ w^Jk/ir. As a result we obtain formulas for the 
hopping 

1 rf2 Vo 



— = / dxwx{x + 1) 



T^ dx^ 2Er 



cos(x) 



w,{x), (27) 



the Raman coupling 

" coii{TTx + (l))\w^{x)\^dx, (/)e{0,7r/2} (28) 



J /**0 / ,/\i / \\2 



Er J Er 
and the interaction energy 



^^ - at"!^ fdx\w,{xr -: at'^hiVo/En). (29) 



Er 

Notice that Wx{x) is now a Wannier function computed for a single, one-dimensional 
problem with period 27r and trapping Vq/ Er. The only other free parameters are 
^o/Er, and the effective one-dimensional scattering length [jj 

at'^^-.^^kat^hiVojERf. (30) 

TT- 

At this point one has to solve a family of one-dimensional problem with the periodic 
potential cos(27ra;) in order to obtain the relative values of J, J' and U. 

5. 2. Experimental parameters 

From the computations sketched in Sect. 15. 1[ it follows that if we choose the Rabi 
coupling of similar strength to the confinement, fio = ^/2, then the ordinary hopping, 
J, and the ladder coupling, J', are also of the same order of magnitude [Fig. [8]. That 
means that by slightly tuning Qq it should be experimentally feasible to cover all 
possible regimes, from J_ — J+ to J_ = 0. 

Now it only remains the question of whether wc can make U/J small or big 
enough to explore the different phases. For a typical experiment with ^^Rb, we have 
a^"^ ~ 5nm and the wavelength of the laser is A ~850nm, so that aj^^, ^ 1.8 x 10~^/^. 
Now, for a strong confinement we have that / ranges from 2 to 8, and comparing 
with the evolution of J and J' in Fig. [5] we see that the value of the interaction can 
indeed change from a superfluid regime, U < \J — J'\, to the Mott insulator regime, 
\J + J'\ -^ U. In other words, if we keep J — J' not too small, the laser intensities will 
not differ much from those used in current experiments [151 139) and we will be able to 
generate all phases. 

While the validity of the tight-binding approximation is well established for the 
ordinary Bose-Hubbard model, wc have introduced a coupling between neighboring 

Ij Note that this is an effective value, unrelated to the scattering theory of cold atoms in ID tubes 
by Olshanii [38j. 
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Figure 8. Parameters of the model, / (solid, left axis), J (dash-dot, right) and 
J' (dashed, right) as a function of the lattice depth, Vo, for a Rabi coupling 
max|n(a:)| = Vb/2, and diagonal interactions, (p = ""/S. The value e (circles, 
right) is the effective coupling to higher bands I I31I I. All parameters are measured 
in units of the recoil energy of the confining lattice. 



sites that might excite the atoms to higher Bloch bands. We have studied this effect 
and computed the foUowing effective strength of the couphng 

2 



£ := 



miix\Q,{no/ER)\ 



AE 



(31) 



which is a function of the energy difference between bands, AE, and the couphng 
between Wannier functions of the lowest and first excited bands 



^{no/ER) := -^ fw^ix)cosix/2 + ^)wi^^ix- 



(32) 



When small, the value e measures how much mixing of the excited wavefunctions there 
is in the ground state and it is a function of fl^. As Fig. [8] shows, for the diagonal 
interactions (0 — tt/2) the coupling is indeed extremely small, e < 10~^, and we can 
be sure that our approximations are valid. Similar results are obtained for = 0. 



5.3. Preparation and observation 

The remaining issues are related to the possibility of preparing and unambiguously 
identifying the ground states of our model ([T]) . The state preparation can be done in a 
two step process. One begins with a ID superfluid and all atoms in the same internal 
state, and slowly increasing the lattice depth. Judging from current experiments, this 
should be a robust way to create a Mott or Tonks gas phase [15] . For U :s> J, J' 
the density excitations will be frozen and the dynamics will be ruled by the effective 
ferromagnetic model p3p . One can can then slowly increase J' to create the ground 
state of (|13l) which will be the ground state of the frustrated model. 

Regarding detection, the three phases of model ([3]) can be identified from time of 
flight images. First of all, for [/ <C J± we expect a quasi-condensate which will produce 
strong interference peaks in either of the atomic species. The width of these peaks may 
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Figure 9. Ideal interference pattern arising from time of flight images of the 
model with diagonal interactions ((S}. We show results for N = L = 32 particles, 
J/U = 0, 0.05 and 0.5, from top to bottom, using J' = J (dash) and J' = 0.9 J 
(solid). See Sect. |531 for details about n and x. 



be slightly modified by temperature |37] but it should remain approximately constant 
for U < J±. For moderately strong interactions, J_ < U < J+, the correlation length 
will decrease substantially as we expect coherences only inside each fragment. Thus 
according to Sect. [33] the interference pattern will be approximately 

|7/.(ip/m)|2 ^ \w{p)\^[l + cos(p)] =: \wip)\^n{p). (33) 

Here w{p) is the Fourier transform of a Wannier wavefunction, n is the ideal 
interference pattern and the position x — tp/m is related to the time, t, the mass 
of the atoms, m, and its momentum in the lattice, p. In practice, for J' ^ J, the 
remaining coherence between fragments can modify this pattern, giving thinner peaks, 
as shown in Fig. [9l 

Another signature of the transition from quasi-condensate to fragmented phase is 
a spontaneous symmetry break that makes the ground state be formed of either c+ or 
c_ particles. The energy gap in this phase is related to the stiffness of the atoms with 
respect to a change in the internal state and it could be probed spectroscopically. 

A third signature which applies to both the fragmented phases and the ground 
state of the Villain model is the existence of magnetic order. More precisely, for strong 
interactions J+ < U, the ground states of the diagonal and Villain frustrated models 
have ferro- and antiferromagnetic order, respectively. This pattern reveals itself in 
high order correlations between the number of atoms of c+ and c_ at different sites: 
{Cj^rCjcr'cl^Cia) oc 1 + eaa' {—iy~^ , where e = 1 or for each model, respectively. 
To observe these correlations one must first rotate the atoms with a 7r/2 pulse, as in 
Eq. ([5]), and afterwards analyze the quantum fluctuations in the time of flight images 
[40] . Alternatively, if the experiment allows for state dependent lattices and photo- 
association, it will be advantageous to use the techniques suggested in [7| to measure 
the number correlations. This second method should indeed provide a much stronger 
signal. 

In current experiments we have to consider two additional sources of imperfection 
that can influence our observations. One of them is the residual harmonic confinement 



Fragmented superfluid due to frustration of cold atoms in optical lattices 18 

which is always present due to the difference on intensity along the lasers that trap the 
atoms. From the diagram of Fig. [5J the harmonic confinement will cause the atomic 
cloud to be made of nested insulating shells with different densities [33], while the 
size of the superfluid regions will substantially decrease and strictly vanish for J = J' . 
This structure can cause noise in the interference patterns, but, as it has been shown 
in recent experiments j41j . it can also be probed using RF knifes and spin-changing 
collisions. 

The second source of imperfections is temperature, which will influence the 
measurements in two ways. First of all, for T ^ 0, the weakly interacting region 
will not be a condensate, but a quasi-condensate, where phase fluctuations are 
energetically cheap but the system remains superfluid. One might argue that these 
phase fluctuations will destroy the interference peaks, but as we have seen in ID 
experiments with optical lattices ^39], this does not seem to happen in practice. 
Moreover, the correlation length of the quasi-condensate will be much larger than 
that of the fragmented phase, a fact that will be evident in the interference pattern. 
It remains the question of whether the fragmented phase itself is robust against 
temperature and we argue that indeed it is. This phase will be observed if the 
temperature of the system is smaller or comparable to the energy gap, A. But as 
shown in Fig. [3]the energy gap can be considerably large, A '^ 0.6 J, sufficiently larger 
than the temperature limitations of the latest experiments [T51 ST] . We expect that 
under these conditions the fragmented phase can be unambiguously identified. 

6. Conclusions 

In this work we have introduced several new ideas. The first one is to use internal 
states of atoms in an optical lattice to simulate additional spatial dimension and to 
implement ladders. While it may be argued that there are previous works which have 
actively employed the internal degrees of freedom of the atom for quantum simulation 
[42] |43] [21 m [6l [7] , it is in this work that they are used to effectively increase the 
dimensionality of the lattice. 

The second idea is to couple the atomic degrees of freedom in a spatially dependent 
way [42], so as to induce a frustrating hopping in the new virtual dimensions. The 
frustrating nature if this assisted tunneling arises from the sign of the Raman coupling 
and it is best appreciated in the hard-core bosons limit, in which the model is 
tantamount to a spin Hamiltonian [Sect. \T/I\ . Additionally, this frustrating hopping 
can also be seen as originating multiple small Mach-Zender interferometers in the 
lattice, such that destructive interference disturb the motion of atoms along the lattice 
[Fig.lB]. 

The most dramatic consequence of these two kinds of tunneling is the breakdown 
of superfluidity. The ground state fragments into a macroscopic number of double- 
well subsystems which lose coherence among them as we increase the strength of 
interactions. This situation is reminiscent of a Bose glass where a disordered energy 
landscape creates islands of superfluid regions with no coherence among them |24| . 
However, in our case the fragmentation is induced by frustration and not by disorder. 

It is possible to qualitatively connect the Physics observed here with that of other 
spin models. On the one hand, the modulation of the hopping leads to an effective 
supperlatice in the internal plus motional degrees of freedom. This superlattice gives 
rise to a larger coherence between neighboring lattice sites, similar to what happens 
in the spin-Peierls effect, where a modulation of a Heisenberg interaction leads to 
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dimerization ^44j . Thus, in contrast to other models where dimerization happens 
spontaneously [22], here it is imposted by the quantum interference of hopping terms. 

However, while this analogy to spin models explains the fragmentation, our system 
goes beyond this simple picture, since we actually have two superlattices which are 
connected by repulsive interactions. The internal degrees of freedom of the atoms 
give rise to rich Physics both in the strongly repulsive regime and in the insulator to 
superfluid transition. In particular, the dynamics of the spin deep in the superfluid 
region deserve tocc be studied in future work. 

The methods shown here can be applied to simulating spin ladders as the ones 
studied in Ref . [7] . Additionally, the internal degrees of freedom of the atoms can be 
used to change the topology of the lattice from square to triangular. Finally, a natural 
extension would be to combine the Raman coupling used here with a two-dimensional 
or even a three-dimensional optical lattice, so as to simulate new topologies and other 
frustrated bilayer magnetic models that may have not been considered so far. 
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